## 
## Night time lights harmonized
##

output_version <- '20201121'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')
st_crs(ntl.zones) <- 4326

# Version 2 Dataset, posted on 08.04.2020, 21:49 by Xuecao Li, Yuyu Zhou, Min zhao, Xia Zhao
# HISTORY
# 17.01.2020 - First online date
# 08.04.2020 - Posted date
# REFERENCES
# https://doi.org/10.3390/rs9060637View

ntl_list = list.files(ntlh_dir,full.names=TRUE,pattern="*.tif")
ntl_list = ntl_list[2:length(ntl_list)]
for (ntl_i in ntl_list){
  r <- raster(ntl_i)
  yyyy <- paste0(substr(names(r),19,22))
  ntl.zones[paste0('ntlharm',yyyy)]<- exact_extract(r,ntl.zones,fun='sum')
  print(paste0("done with ",yyyy))
}

out.df <- as.data.frame(st_drop_geometry(ntl.zones))
save.dta13(out.df,paste0(wkdir,'/proc_data/NTLHARM_ts',output_version,'.dta'))